package flexibilityDesign;

import ilog.concert.IloConversion;
import ilog.concert.IloLinearNumExpr;
import ilog.concert.IloNumVar;
import ilog.concert.IloNumVarType;
import ilog.concert.IloRange;
import ilog.cplex.IloCplex;

import java.io.BufferedWriter;
import java.io.FileWriter;
import java.io.IOException;
import java.io.Writer;
import java.math.BigDecimal;
import java.util.ArrayList;
import java.util.HashMap;
import java.util.List;
import java.util.Map;
import java.util.Stack;

import concurrent.Executer;
import support.Aggregate;
import support.AlgorithmSupport;
import util.Lang;
import util.NullOutputStream;

public class StructureDesignHybridCut {
	private int product_num;// product number
	private int plant_num;// plant number
	private int simulationCount;// sampling number

	private int groups = 80;// groups count

	// enhancement flag
	private boolean parallelComput = false;
	private boolean eliminateNonImprove = false;
	private boolean trustRegion = false;
	private boolean approMasterSolve = false;
	private boolean initSolve = false;
	private boolean strengtheningCut = false;

	private int simulation_index = 0;

	public void setSimulationIndex(int index) {
		this.simulation_index = index;
	}

	private Map<Integer, List<Double[]>> demands;// clustered demands
	private double[] capacities;// capacities of plants
	private double[][] link_c;// link cost
	private double[][] profit_m;// profit magin

	public StructureDesignHybridCut(int n, int m, int simulationC, int groups) {
		this.product_num = n;
		this.plant_num = m;
		this.simulationCount = simulationC;
		this.groups = groups;
	}

	public void setEnhances(boolean paral, boolean eliminateNonImprove, boolean trustRegion, boolean approMaster,
			boolean initSolve, boolean strengthenCut) {
		this.parallelComput = paral;
		this.eliminateNonImprove = eliminateNonImprove;
		this.trustRegion = trustRegion;
		this.approMasterSolve = approMaster;
		this.initSolve = initSolve;
		this.strengtheningCut = strengthenCut;
	}

	private void initData() {
		System.out.println("Init data...");
		this.capacities = Support.readCapacities(product_num, plant_num);

		// double[][] rawDemands = Support.readDemands(simulationCount,
		// product_num, simulation_index);
		double[][] rawDemands = Support.readDemands(simulationCount, product_num);

		Aggregate aggregate = new Aggregate();
		// this.demands = aggregate.kmeans(rawDemands, this.groups);
		// this.demands = aggregate.roundRobin(rawDemands, groups);
		this.demands = aggregate.kmeans_matlab(rawDemands, this.groups);

		// this.demands = Support.readDemands_clustered(simulationCount, n);

		this.link_c = Support.readLinkCost(product_num, plant_num);
		this.profit_m = Support.readProfitMargin(product_num, plant_num);

		System.out.println("Finish init data");
	}

	public int[][] Lshape() {
		this.initData();

		List<Double[]> bounds = new ArrayList<>();

		// probability of total
		BigDecimal pro_t = BigDecimal.ONE.divide(BigDecimal.valueOf(simulationCount), 5, BigDecimal.ROUND_HALF_EVEN);

		int G = this.demands.size();
		// probability of groups
		BigDecimal pro_g = BigDecimal.ONE.divide(BigDecimal.valueOf(G), 5, BigDecimal.ROUND_HALF_EVEN);

		double LB = Double.NEGATIVE_INFINITY;
		double UB = Double.MAX_VALUE;
		double gap = Double.MAX_VALUE;

		Stack<Double> historyLB = new Stack<>();

		boolean stable = false;// drop trust region once stable
		double preGap = Double.MAX_VALUE;

		long timeLimite = AlgorithmSupport.TERMINATE_TIME_LIMITE;

		try {
			IloCplex master = new IloCplex();
			master.setOut(new NullOutputStream());
			IloNumVar[][] Y = new IloNumVar[product_num][plant_num];
			int[][] YValue = new int[product_num][plant_num];
			IloNumVar[] thetas = master.numVarArray(G, -1000000, 1000000);
			double[] thetasValue = new double[G];
			for (int i = 0; i < product_num; i++) {
				for (int j = 0; j < plant_num; j++) {
					Y[i][j] = master.boolVar();
				}
			}

			// @enhance warm start up from small sample solutions
			if (this.initSolve) {
				DesignMIP dmip = new DesignMIP(product_num, plant_num, 50, 5);
				dmip.turnOffWriteToFile();
				int[][] warmY = dmip.MIPSolve();
				YValue = Lang.clone(warmY);
			} else {
				for (int i = 0; i < product_num; i++)
					for (int j = 0; j < plant_num; j++)
						YValue[i][j] = 1;
			}
			// end enhance

			for (int g = 0; g < G; g++)
				thetasValue[g] = -100000;// init theta value

			double[][] corePoint = null;

			// all set to be 1
			corePoint = new double[product_num][plant_num];
			for (int i = 0; i < product_num; i++)
				for (int j = 0; j < plant_num; j++)
					corePoint[i][j] = 1;

			IloLinearNumExpr masterObj = master.linearNumExpr();
			for (int i = 0; i < product_num; i++) {
				for (int j = 0; j < plant_num; j++) {
					masterObj.addTerm(link_c[i][j], Y[i][j]);
				}
			}
			for (int g = 0; g < G; g++) {
				// masterObj.addTerm(pro_g.doubleValue(), thetas[g]);
				masterObj.addTerm(1.0, thetas[g]);
			}
			master.addMinimize(masterObj);

			int[][] currentOptSolution = Lang.clone(YValue);
			double currentOptVal = 0.0;
			int delta = product_num * plant_num / 2;
			int samegapCount = 0;
			double currentGap = Double.MAX_VALUE;

			long start = System.currentTimeMillis();
			int iter = 1;

			while (true) {

				if (iter > 1) {// solve master problem
					// @enhance: dynamically adjust tolerance
					if (this.approMasterSolve) {
						if (samegapCount >= 2) {
							double gapT = master.getParam(IloCplex.DoubleParam.EpGap);
							gapT = gapT / 2;
							if (gapT < 0.001)
								gapT = 0.001;
							master.setParam(IloCplex.DoubleParam.EpGap, gapT);
						} else {
							if (gap > 0.05)
								master.setParam(IloCplex.DoubleParam.EpGap, 0.05);
							// else if (gap > 0.03)
							// master.setParam(IloCplex.DoubleParam.EpGap,
							// 0.03);
							// else if (gap > 0.02)
							// master.setParam(IloCplex.DoubleParam.EpGap,
							// 0.08);
							else if (gap > 0.01)
								master.setParam(IloCplex.DoubleParam.EpGap, 0.01);
							// else if (gap > 0.008)
							// master.setParam(IloCplex.DoubleParam.EpGap,
							// 0.008);
							// else if (gap > 0.007)
							// master.setParam(IloCplex.DoubleParam.EpGap,
							// 0.007);
							else if (gap > 0.005)
								master.setParam(IloCplex.DoubleParam.EpGap, 0.005);
							else if (gap > 0.003)
								master.setParam(IloCplex.DoubleParam.EpGap, 0.003);
							// else if (gap > 0.002)
							// master.setParam(IloCplex.DoubleParam.EpGap,
							// 0.005);
							else
								master.setParam(IloCplex.DoubleParam.EpGap, 0.001);
						}
					}
					// end enhance

					// @enhance: trust region
					IloRange trustRegion = null;
					if (this.trustRegion && !stable) {
						IloLinearNumExpr trustExpr = master.linearNumExpr();
						int onecount = 0;
						for (int i = 0; i < product_num; i++) {
							for (int j = 0; j < plant_num; j++) {
								if (currentOptSolution[i][j] == 1) {
									onecount++;
									trustExpr.addTerm(-1, Y[i][j]);
								} else
									trustExpr.addTerm(1, Y[i][j]);
							}
						}
						trustRegion = master.addLe(trustExpr, delta - onecount);
					}
					// end enhance

					// master.exportModel("D:\\master_" + iter + ".lp");

					for (int i = 0; i < product_num; i++) {
						double[] init = new double[plant_num];
						for (int j = 0; j < plant_num; j++)
							init[j] = currentOptSolution[i][j];
						master.addMIPStart(Y[i], init);
					}

					master.setParam(IloCplex.DoubleParam.TiLim, timeLimite/1000);

					long masterStart = System.currentTimeMillis();
					if (master.solve()) {
						long masterEnd = System.currentTimeMillis();
						System.out.println("master time = " + (masterEnd - masterStart) + "; total time="
								+ (masterEnd - start));

						double currentMasterObjV = master.getObjValue();
						if (currentMasterObjV > LB) {
							historyLB.add(LB);
							LB = currentMasterObjV;
						}

						// dynamic changing delta
						if (!stable && this.trustRegion) {
							if (currentMasterObjV > LB)
								delta = (int) Math.ceil(delta * 1.3);
							else if (iter > 2 && gap > 0)
								delta = (int) Math.ceil(delta * 0.6);
							if (delta < 1)
								delta = 1;
						}
						System.out.println("delta=" + delta);

						for (int i = 0; i < product_num; i++)
							for (int j = 0; j < plant_num; j++)
								YValue[i][j] = (int) master.getValue(Y[i][j]);
						thetasValue = master.getValues(thetas);
					}

					// @enhance strengthening cuts
					long strengtheningStart = System.currentTimeMillis();
					if (this.strengtheningCut) {
						corePoint = new double[product_num][plant_num];
						IloConversion[][] conversions = new IloConversion[product_num][plant_num];
						for (int i = 0; i < product_num; i++) {
							for (int j = 0; j < plant_num; j++) {
								conversions[i][j] = master.conversion(Y[i][j], IloNumVarType.Float);
								master.add(conversions[i][j]);
							}
						}
						if (master.solve()) {
							for (int i = 0; i < product_num; i++)
								corePoint[i] = master.getValues(Y[i]);
						}
						for (int i = 0; i < product_num; i++) {
							for (int j = 0; j < plant_num; j++) {
								master.remove(conversions[i][j]);
							}
						}
					}
					long strengtheningEnd = System.currentTimeMillis();
					System.out.println("@relaxed solve time = " + (strengtheningEnd - strengtheningStart));

					if (trustRegion != null)
						master.delete(trustRegion);
				}

				double estimatedValue = 0.0;// extensive form value
				double subObjValue = 0.0;
				for (int i = 0; i < product_num; i++) {
					for (int j = 0; j < plant_num; j++) {
						if (YValue[i][j] == 1)
							estimatedValue += link_c[i][j];
					}
				}

				boolean hasMoreCut = false;
				long subStart = System.currentTimeMillis();
				for (int g = 0; g < G; g++) {// handle each group
					List<Double[]> demand_g = demands.get(g);
					int SG = demand_g.size();// sub group size

					BigDecimal et = BigDecimal.ZERO;
					BigDecimal[][] Et = new BigDecimal[product_num][plant_num];
					for (int i = 0; i < product_num; i++) {
						for (int j = 0; j < plant_num; j++)
							Et[i][j] = BigDecimal.ZERO;
					}

					BigDecimal pro_s = BigDecimal.ONE.divide(BigDecimal.valueOf(SG), 5, BigDecimal.ROUND_HALF_EVEN);

					BigDecimal wv = BigDecimal.ZERO;

					if (parallelComput) {
						// @enhance parallel computing
						Executer executer = new Executer(10, true);
						for (int k = 0; k < SG; k++) {
							// @enhance strengthen cut
							if (this.strengtheningCut)
								executer.fork(new DualSolver(product_num, plant_num, Lang.convertDouble2double(demand_g
										.get(k)), YValue, corePoint, profit_m, capacities, k));
							else
								executer.fork(new DualSolver(product_num, plant_num, Lang.convertDouble2double(demand_g
										.get(k)), YValue, profit_m, capacities, k));
						}
						List<DualResult> results = executer.join();
						executer.shutdown();
						for (DualResult dr : results) {
							if (dr != null) {
								subObjValue += dr.getObjValue();
								et = et.add(dr.getEtm());
								for (int i = 0; i < product_num; i++)
									for (int j = 0; j < plant_num; j++)
										Et[i][j] = Et[i][j].add(dr.getEt()[i][j]);
							}
						}
						// enhance parallel computing end
					} else {
						for (int k = 0; k < SG; k++) {// handle each scenario
							IloCplex dual = new IloCplex();
							dual.setOut(new NullOutputStream());
							IloNumVar[] phi = dual.numVarArray(plant_num, Double.NEGATIVE_INFINITY, 0.0);
							IloNumVar[] sigma = dual.numVarArray(product_num, Double.NEGATIVE_INFINITY, 0.0);
							IloNumVar[][] gamma = new IloNumVar[product_num][plant_num];
							for (int i = 0; i < product_num; i++)
								for (int j = 0; j < plant_num; j++)
									gamma[i][j] = dual.numVar(Double.NEGATIVE_INFINITY, 0.0);

							IloLinearNumExpr dualObj = dual.linearNumExpr();
							for (int j = 0; j < plant_num; j++)
								dualObj.addTerm(capacities[j], phi[j]);
							for (int i = 0; i < product_num; i++)
								dualObj.addTerm(demand_g.get(k)[i], sigma[i]);
							for (int i = 0; i < product_num; i++) {
								for (int j = 0; j < plant_num; j++) {
									dualObj.addTerm(YValue[i][j] * capacities[j], gamma[i][j]);

									// constraint
									dual.addLe(dual.sum(phi[j], sigma[i], gamma[i][j]), -profit_m[i][j]);
								}
							}
							dual.addMaximize(dualObj);

							if (dual.solve()) {
								double objv = dual.getObjValue();

								subObjValue += objv;
								for (int i = 0; i < product_num; i++) {
									et = et.add(BigDecimal.valueOf(dual.getValue(sigma[i]) * demand_g.get(k)[i]));
								}
								for (int j = 0; j < plant_num; j++) {
									et = et.add(BigDecimal.valueOf(dual.getValue(phi[j]) * capacities[j]));
								}
								for (int i = 0; i < product_num; i++) {
									for (int j = 0; j < plant_num; j++) {
										Et[i][j] = Et[i][j].add(BigDecimal.valueOf(dual.getValue(gamma[i][j])
												* (-capacities[j])));
									}
								}
							}
							dual.end();
						}
					}// end for one group
						// et = et.multiply(pro_s);
					et = et.multiply(pro_t);
					for (int i = 0; i < product_num; i++) {
						for (int j = 0; j < plant_num; j++) {
							// Et[i][j] = Et[i][j].multiply(pro_s);
							Et[i][j] = Et[i][j].multiply(pro_t);
						}
					}
					BigDecimal EtY = BigDecimal.ZERO;
					for (int i = 0; i < product_num; i++) {
						for (int j = 0; j < plant_num; j++) {
							if (YValue[i][j] == 1)
								EtY.add(Et[i][j]);
						}
					}
					wv = et.subtract(EtY);

					if (wv.doubleValue() > thetasValue[g]) {// generate cut
						hasMoreCut = true;

						IloLinearNumExpr cut = master.linearNumExpr();
						for (int i = 0; i < product_num; i++) {
							for (int j = 0; j < plant_num; j++) {
								cut.addTerm(Et[i][j].doubleValue(), Y[i][j]);
							}
						}
						cut.addTerm(1.0, thetas[g]);
						// add optimality cut
						master.addGe(cut, et.doubleValue());
					}
				}// end for all groups
				long subEnd = System.currentTimeMillis();
				System.out.println("@subproblem time=" + (subEnd - subStart));
				estimatedValue += BigDecimal.valueOf(subObjValue).multiply(pro_t).doubleValue();
				if (estimatedValue < UB)
					UB = estimatedValue;

				// @enhance eliminate current solution
				if (estimatedValue < currentOptVal) {
					currentOptVal = estimatedValue;
					currentOptSolution = Lang.clone(YValue);
				} else if (this.eliminateNonImprove) {
					if (iter > 3) {
						IloLinearNumExpr cutExp = master.linearNumExpr();
						int countOne = 0;
						for (int i = 0; i < product_num; i++) {
							for (int j = 0; j < plant_num; j++) {
								if (YValue[i][j] == 1) {
									countOne++;
									cutExp.addTerm(-1, Y[i][j]);
								} else
									cutExp.addTerm(1, Y[i][j]);
							}
						}
						master.addGe(cutExp, 1 - countOne);
					}
				}
				// end enhance

				preGap = gap;

				if (iter > 1) {
					gap = BigDecimal.valueOf(UB - LB).divide(BigDecimal.valueOf(-LB), 5, BigDecimal.ROUND_HALF_EVEN)
							.doubleValue();
				} else {
					gap = UB - LB;
				}

				if (gap == currentGap) {
					samegapCount++;
				} else {
					samegapCount = 0;
					currentGap = gap;
				}

				System.out.println("stable gap = " + (preGap - gap));

				if (preGap - gap < 0.004 && samegapCount < 2)
					stable = true;
				else
					stable = false;

				long end = System.currentTimeMillis();
				long time = end - start;
				System.out.println("iter=" + iter + "\t UB=" + UB + "\t LB=" + LB + "\t gap=" + gap);

				Double[] bound = new Double[2];
				bound[0] = LB;
				bound[1] = UB;
				bounds.add(bound);

				timeLimite = AlgorithmSupport.TERMINATE_TIME_LIMITE - time;

				System.out.println("timelimit left: " + timeLimite + "; time=" + time);

				if ((gap < 0.001 && gap >= 0.000) || !hasMoreCut || time > AlgorithmSupport.TERMINATE_TIME_LIMITE) {
					System.out.println((gap < 0.001 && gap >= 0.000) + "," + (!hasMoreCut) + ","
							+ AlgorithmSupport.TERMINATE_TIME_LIMITE);

					System.out.println("hasMorecut=" + hasMoreCut);
					StringBuilder builder = new StringBuilder("\r\n n,m,M=" + product_num + "," + plant_num + ","
							+ simulationCount + "\r\n" + "HybridCut \t" + "group=" + G + "; index=" + simulation_index
							+ "\r\n");

					if (this.parallelComput)
						builder.append("parallel \t");
					if (this.eliminateNonImprove)
						builder.append("eliminate non improve \t");
					if (this.trustRegion)
						builder.append("trust region \t");
					if (this.approMasterSolve)
						builder.append("approximate master solve \t");
					if (this.initSolve)
						builder.append("init solver \t");
					if (this.strengtheningCut)
						builder.append("strengthen cuts \t");

					builder.append("\r\n");
					for (int i = 0; i < product_num; i++) {
						for (int j = 0; j < plant_num; j++) {
							if (YValue[i][j] == 1)
								System.out.print("(" + i + "," + j + ")\t");
						}
						System.out.println();
					}
					double value = BigDecimal.valueOf(LB + UB)
							.divide(BigDecimal.valueOf(2), 5, BigDecimal.ROUND_HALF_EVEN).doubleValue();

					System.out.println("time consumed: " + time);
					builder.append("time=" + time + "\t" + "gap=" + gap + "\t value=" + value + "\r\n");
					for (Double[] bo : bounds)
						builder.append(bo[0] + "," + bo[1] + "\r\n");

					// try {
					// Writer writer = new BufferedWriter(new
					// FileWriter(Lang.projectRoot()
					// + "\\data\\flexibilityDesign\\result\\result.txt",
					// true));
					// writer.write(builder.toString());
					// writer.close();
					// } catch (IOException e) {
					// e.printStackTrace();
					// }
					break;
				} else {
					iter++;
					if (gap < 0.0) {
						delta = 2 * delta;
						if (delta > product_num * plant_num)
							delta = product_num * plant_num;
						if (historyLB.size() > 0)
							LB = historyLB.pop();
					}
				}
			}
			master.end();

			return YValue;
		} catch (Exception e) {
			e.printStackTrace();
			return null;
		}
	}
}
